﻿#include <math.h>
#include <malloc.h>

/**
 * @brief cagaus 全选主元高斯消去法
 * @param a n*n 存放方程组的系数矩阵，返回时将被破坏
 * @param b 常数向量
 * @param n 方程组的阶数
 * @param x 返回方程组的解向量
 * @param perm 辅助数组，n 个元素。函数利用这个数组来记录选主元过程中的位置交换。
 * @return 返回0表示原方程组的系数矩阵奇异
 */
int cagaus(double a[], double b[], int n, double x[], int perm[])
{
    int l, k, i, j, is, p, q;
    double d, t;
    l = 1;
    for(k = 0; k <= n - 2; k++)
    {
        d = 0.0;
        for(i = k; i <= n - 1; i++)
        {
            for(j = k; j <= n - 1; j++)
            {
                t = fabs(a[i * n + j]);
                if(t > d)
                {
                    d = t;
                    perm[k] = j;
                    is = i;
                }
            }
        }
        if(d + 1.0 == 1.0)
        {
            l = 0;
        }
        else
        {
            if(perm[k] != k)
            {
                for(i = 0; i <= n - 1; i++)
                {
                    p = i * n + k;
                    q = i * n + perm[k];
                    t = a[p];
                    a[p] = a[q];
                    a[q] = t;
                }
            }
            if(is != k)
            {
                for(j = k; j <= n - 1; j++)
                {
                    p = k * n + j;
                    q = is * n + j;
                    t = a[p];
                    a[p] = a[q];
                    a[q] = t;
                }
                t = b[k];
                b[k] = b[is];
                b[is] = t;
            }
        }
        if(l == 0)
        {
            return(0);
        }
        d = a[k * n + k];
        for(j = k + 1; j <= n - 1; j++)
        {
            p = k * n + j;
            a[p] = a[p] / d;
        }
        b[k] = b[k] / d;
        for(i = k + 1; i <= n - 1; i++)
        {
            for(j = k + 1; j <= n - 1; j++)
            {
                p = i * n + j;
                a[p] = a[p] - a[i * n + k] * a[k * n + j];
            }
            b[i] = b[i] - a[i * n + k] * b[k];
        }
    }
    d = a[(n - 1) * n + n - 1];
    if(fabs(d) + 1.0 == 1.0)
    {
        return(0);
    }
    x[n - 1] = b[n - 1] / d;
    for(i = n - 2; i >= 0; i--)
    {
        t = 0.0;
        for(j = i + 1; j <= n - 1; j++)
        {
            t = t + a[i * n + j] * x[j];
        }
        x[i] = b[i] - t;
    }
    perm[n - 1] = n - 1;
    for(k = n - 1; k >= 0; k--)
    {
        if(perm[k] != k)
        {
            t = x[k];
            x[k] = x[perm[k]];
            x[perm[k]] = t;
        }
    }

    return 1;
}

/**
  * @brief cbcgaus 复系数方程组的全选主元高斯消去法
  * @param ar n * n存放复系数矩阵的实部，要被破坏
  * @param ai n * n存放复系数矩阵的虚部，要被破坏
  * @param n 方程组的阶数
  * @param br 存放方程组右端复常数向量的实部，返回时存放解向量的实部
  * @param bi 存放方程组右端复常数向量的虚部，返回时存放解向量的虚部
  * @param perm 辅助数组，n 个元素。函数利用这个数组来记录选主元过程中的位置交换。
  * @return  返回值0表示系数矩阵奇异
  */
int cbcgaus(double ar[], double ai[], int n, double br[], double bi[], int perm[])
{
    int l, k, i, j, is, u, v;
    double p, q, s, d;
    for(k = 0; k <= n - 2; k ++)
    {
        d = 0.0;
        for(i = k; i <= n - 1; i ++)
        {
            for(j = k; j <= n - 1; j ++)
            {
                u = i * n + j;
                p = ar[u] * ar[u] + ai[u] * ai[u];
                if(p > d)
                {
                    d = p;
                    perm[k] = j;
                    is = i;
                }
            }
        }
        if(d + 1.0 == 1.0)
        {
            return(0);
        }
        if(is != k)
        {
            for(j = k; j <= n - 1; j ++)
            {
                u = k * n + j;
                v = is * n + j;
                p = ar[u];
                ar[u] = ar[v];
                ar[v] = p;
                p = ai[u];
                ai[u] = ai[v];
                ai[v] = p;
            }
            p = br[k];
            br[k] = br[is];
            br[is] = p;
            p = bi[k];
            bi[k] = bi[is];
            bi[is] = p;
        }
        if(perm[k] != k)
        {
            for(i = 0; i <= n - 1; i ++)
            {
                u = i * n + k;
                v = i * n + perm[k];
                p = ar[u];
                ar[u] = ar[v];
                ar[v] = p;
                p = ai[u];
                ai[u] = ai[v];
                ai[v] = p;
            }
        }
        v = k * n + k;
        for(j = k + 1; j <= n - 1; j ++)
        {
            u = k * n + j;
            p = ar[u] * ar[v];
            q = - ai[u] * ai[v];
            s = (ar[v] - ai[v]) * (ar[u] + ai[u]);
            ar[u] = (p - q) / d;
            ai[u] = (s - p - q) / d;
        }
        p = br[k] * ar[v];
        q = - bi[k] * ai[v];
        s = (ar[v] - ai[v]) * (br[k] + bi[k]);
        br[k] = (p - q) / d;
        bi[k] = (s - p - q) / d;
        for(i = k + 1; i <= n - 1; i ++)
        {
            u = i * n + k;
            for(j = k + 1; j <= n - 1; j ++)
            {
                v = k * n + j;
                l = i * n + j;
                p = ar[u] * ar[v];
                q = ai[u] * ai[v];
                s = (ar[u] + ai[u]) * (ar[v] + ai[v]);
                ar[l] = ar[l] - p + q;
                ai[l] = ai[l] - s + p + q;
            }
            p = ar[u] * br[k];
            q = ai[u] * bi[k];
            s = (ar[u] + ai[u]) * (br[k] + bi[k]);
            br[i] = br[i] - p + q;
            bi[i] = bi[i] - s + p + q;
        }
    }
    u = (n - 1) * n + n - 1;
    d = ar[u] * ar[u] + ai[u] * ai[u];
    if(d + 1.0 == 1.0)
    {
        return(0);
    }
    p = ar[u] * br[n - 1];
    q = - ai[u] * bi[n - 1];
    s = (ar[u] - ai[u]) * (br[n - 1] + bi[n - 1]);
    br[n - 1] = (p - q) / d;
    bi[n - 1] = (s - p - q) / d;
    for(i = n - 2; i >= 0; i --)
    {
        for(j = i + 1; j <= n - 1; j ++)
        {
            u = i * n + j;
            p = ar[u] * br[j];
            q = ai[u] * bi[j];
            s = (ar[u] + ai[u]) * (br[j] + bi[j]);
            br[i] = br[i] - p + q;
            bi[i] = bi[i] - s + p + q;
        }
    }
    perm[n - 1] = n - 1;
    for(k = n - 1; k >= 0; k --)
    {
        if(perm[k] != k)
        {
            p = br[k];
            br[k] = br[perm[k]];
            br[perm[k]] = p;
            p = bi[k];
            bi[k] = bi[perm[k]];
            bi[perm[k]] = p;
        }
    }
    return(1);
}

/**
 * @brief ccgj 全选主元高斯-约当消去法
 * @param a  n*n 存放方程组的系数矩阵，返回时将被破坏
 * @param b  n*m常数向量,每列为一组，返回时存放m组解向量
 * @param n  方程组的阶数
 * @param m  方程组右端常数向量的个数
 * @param perm 辅助数组，n 个元素。函数利用这个数组来记录选主元过程中的位置交换。
 * @return 返回0表示原方程组的系数矩阵奇异
 */
int ccgj(double a[], double b[], int n, int m, int perm[])
{
    int l, k, i, j, is, p, q;
    double d, t;
    l = 1;
    for(k = 0; k <= n - 1; k++)
    {
        d = 0.0;
        for(i = k; i <= n - 1; i++)
        {
            for(j = k; j <= n - 1; j++)
            {
                t = fabs(a[i * n + j]);
                if(t > d)
                {
                    d = t;
                    perm[k] = j;
                    is = i;
                }
            }
        }
        if(d + 1.0 == 1.0)
        {
            l = 0;
        }
        else
        {
            if(perm[k] != k)
            {
                for(i = 0; i <= n - 1; i++)
                {
                    p = i * n + k;
                    q = i * n + perm[k];
                    t = a[p];
                    a[p] = a[q];
                    a[q] = t;
                }
            }
            if(is != k)
            {
                for(j = k; j <= n - 1; j++)
                {
                    p = k * n + j;
                    q = is * n + j;
                    t = a[p];
                    a[p] = a[q];
                    a[q] = t;
                }
                for(j = 0; j <= m - 1; j++)
                {
                    p = k * m + j;
                    q = is * m + j;
                    t = b[p];
                    b[p] = b[q];
                    b[q] = t;
                }
            }
        }
        if(l == 0)
        {
            return(0);
        }
        d = a[k * n + k];
        for(j = k + 1; j <= n - 1; j++)
        {
            p = k * n + j;
            a[p] = a[p] / d;
        }
        for(j = 0; j <= m - 1; j++)
        {
            p = k * m + j;
            b[p] = b[p] / d;
        }
        for(j = k + 1; j <= n - 1; j++)
        {
            for(i = 0; i <= n - 1; i++)
            {
                p = i * n + j;
                if(i != k)
                {
                    a[p] = a[p] - a[i * n + k] * a[k * n + j];
                }
            }
        }
        for(j = 0; j <= m - 1; j++)
        {
            for(i = 0; i <= n - 1; i++)
            {
                p = i * m + j;
                if(i != k)
                {
                    b[p] = b[p] - a[i * n + k] * b[k * m + j];
                }
            }
        }
    }
    for(k = n - 1; k >= 0; k--)
    {
        if(perm[k] != k)
        {
            for(j = 0; j <= m - 1; j++)
            {
                p = k * m + j;
                q = perm[k] * m + j;
                t = b[p];
                b[p] = b[q];
                b[q] = t;
            }
        }
    }
    return(1);
}

/**
 * @brief cdcgj 复系数全选主元高斯-约当消去法
 * @param ar n*n 存放方程组的系数矩阵的实部，返回时将被破坏
 * @param ai n*n 存放方程组的系数矩阵的虚部，返回时将被破坏
 * @param br n*m常数向量的实部,每列为一组，返回时存放m组解向量的实部
 * @param bi n*m常数向量的虚部,每列为一组，返回时存放m组解向量的虚部
 * @param n 方程组的阶数
 * @param m 方程组右端常数向量的个数
 * @param perm 辅助数组，n 个元素。函数利用这个数组来记录选主元过程中的位置交换。
 * @return 返回0表示原方程组的系数矩阵奇异
 */
int cdcgj(double ar[], double ai[], double br[], double bi[], int n, int m, int perm[])
{
    int l, k, i, j, is, u, v;
    double p, q, s, d;
    for(k = 0; k <= n - 1; k++)
    {
        d = 0.0;
        for(i = k; i <= n - 1; i++)
        {
            for(j = k; j <= n - 1; j++)
            {
                u = i * n + j;
                p = ar[u] * ar[u] + ai[u] * ai[u];
                if(p > d)
                {
                    d = p;
                    perm[k] = j;
                    is = i;
                }
            }
        }
        if(d + 1.0 == 1.0)
        {
            return(0);
        }
        if(is != k)
        {
            for(j = k; j <= n - 1; j++)
            {
                u = k * n + j;
                v = is * n + j;
                p = ar[u];
                ar[u] = ar[v];
                ar[v] = p;
                p = ai[u];
                ai[u] = ai[v];
                ai[v] = p;
            }
            for(j = 0; j <= m - 1; j++)
            {
                u = k * m + j;
                v = is * m + j;
                p = br[u];
                br[u] = br[v];
                br[v] = p;
                p = bi[u];
                bi[u] = bi[v];
                bi[v] = p;
            }
        }
        if(perm[k] != k)
        {
            for(i = 0; i <= n - 1; i++)
            {
                u = i * n + k;
                v = i * n + perm[k];
                p = ar[u];
                ar[u] = ar[v];
                ar[v] = p;
                p = ai[u];
                ai[u] = ai[v];
                ai[v] = p;
            }
        }
        v = k * n + k;
        for(j = k + 1; j <= n - 1; j++)
        {
            u = k * n + j;
            p = ar[u] * ar[v];
            q = -ai[u] * ai[v];
            s = (ar[v] - ai[v]) * (ar[u] + ai[u]);
            ar[u] = (p - q) / d;
            ai[u] = (s - p - q) / d;
        }
        for(j = 0; j <= m - 1; j++)
        {
            u = k * m + j;
            p = br[u] * ar[v];
            q = -bi[u] * ai[v];
            s = (ar[v] - ai[v]) * (br[u] + bi[u]);
            br[u] = (p - q) / d;
            bi[u] = (s - p - q) / d;
        }
        for(i = 0; i <= n - 1; i++)
        {
            if(i != k)
            {
                u = i * n + k;
                for(j = k + 1; j <= n - 1; j++)
                {
                    v = k * n + j;
                    l = i * n + j;
                    p = ar[u] * ar[v];
                    q = ai[u] * ai[v];
                    s = (ar[u] + ai[u]) * (ar[v] + ai[v]);
                    ar[l] = ar[l] - p + q;
                    ai[l] = ai[l] - s + p + q;
                }
                for(j = 0; j <= m - 1; j++)
                {
                    l = i * m + j;
                    v = k * m + j;
                    p = ar[u] * br[v];
                    q = ai[u] * bi[v];
                    s = (ar[u] + ai[u]) * (br[v] + bi[v]);
                    br[l] = br[l] - p + q;
                    bi[l] = bi[l] - s + p + q;
                }
            }
        }
    }
    for(k = n - 1; k >= 0; k--)
    {
        if(perm[k] != k)
        {
            for(j = 0; j <= m - 1; j++)
            {
                u = k * m + j;
                v = perm[k] * m + j;
                p = br[u];
                br[u] = br[v];
                br[v] = p;
                p = bi[u];
                bi[u] = bi[v];
                bi[v] = p;
            }
        }
    }
    return(1);
}

/**
 * @brief cetrd  求解三对角线方程组的追赶法
 * @param b 存放三对角线矩阵中三对角线上的元素,
 *          存放顺序为a00, a01, a10, a11, a12, a21, a22, a23 ,..., an-1, n-2, an-1, n-1, 就是按行依次存储
 * @param n 方程组的阶数
 * @param m b 的长度，应为 m = 3n - 2, 本函数中用来做检验
 * @param d 存放方程组右端常数向量，返回时存放解向量
 * @return 1 表示成功，其他表示失败
 */
int cetrd(double b[], int n, int m, double d[])
{
    int k, j;
    double s;
    if(m != (3 * n - 2))
    {
        return(-2);
    }
    for(k = 0; k <= n - 2; k++)
    {
        j = 3 * k;
        s = b[j];
        if(fabs(s) + 1.0 == 1.0)
        {
            return(0);
        }
        b[j + 1] = b[j + 1] / s;
        d[k] = d[k] / s;
        b[j + 3] = b[j + 3] - b[j + 2] * b[j + 1];
        d[k + 1] = d[k + 1] - b[j + 2] * d[k];
    }
    s = b[3 * n - 3];
    if(fabs(s) + 1.0 == 1.0)
    {
        return(0);
    }
    d[n - 1] = d[n - 1] / s;
    for(k = n - 2; k >= 0; k--)
    {
        d[k] = d[k] - b[3 * k + 1] * d[k + 1];
    }
    return(2);
}

/**
 * @brief cfband 一般带型方程组的求解
 * @param b
 * @param d 存放方程组右端常数向量，返回时存放解向量
 * @param n 方程组的阶数
 * @param l 系数矩阵的半带宽
 * @param il 系数矩阵的带宽，应为il=2l+1
 * @param m 方程组右端的常数
 * @return
 */
int cfband(double b[], double d[], int n, int l, int il, int m)
{
    int ls, k, i, j, is, u, v;
    double p, t;
    if(il != (2 * l + 1))
    {
        return(-2);
    }
    ls = l;
    for(k = 0; k <= n - 2; k++)
    {
        p = 0.0;
        for(i = k; i <= ls; i++)
        {
            t = fabs(b[i * il]);
            if(t > p)
            {
                p = t;
                is = i;
            }
        }
        if(p + 1.0 == 1.0)
        {
            return(0);
        }
        for(j = 0; j <= m - 1; j++)
        {
            u = k * m + j;
            v = is * m + j;
            t = d[u];
            d[u] = d[v];
            d[v] = t;
        }
        for(j = 0; j <= il - 1; j++)
        {
            u = k * il + j;
            v = is * il + j;
            t = b[u];
            b[u] = b[v];
            b[v] = t;
        }
        for(j = 0; j <= m - 1; j++)
        {
            u = k * m + j;
            d[u] = d[u] / b[k * il];
        }
        for(j = 1; j <= il - 1; j++)
        {
            u = k * il + j;
            b[u] = b[u] / b[k * il];
        }
        for(i = k + 1; i <= ls; i++)
        {
            t = b[i * il];
            for(j = 0; j <= m - 1; j++)
            {
                u = i * m + j;
                v = k * m + j;
                d[u] = d[u] - t * d[v];
            }
            for(j = 1; j <= il - 1; j++)
            {
                u = i * il + j;
                v = k * il + j;
                b[u - 1] = b[u] - t * b[v];
            }
            u = i * il + il - 1;
            b[u] = 0.0;
        }
        if(ls != (n - 1))
        {
            ls = ls + 1;
        }
    }
    p = b[(n - 1) * il];
    if(fabs(p) + 1.0 == 1.0)
    {
        return(0);
    }
    for(j = 0; j <= m - 1; j++)
    {
        u = (n - 1) * m + j;
        d[u] = d[u] / p;
    }
    ls = 1;
    for(i = n - 2; i >= 0; i--)
    {
        for(k = 0; k <= m - 1; k++)
        {
            u = i * m + k;
            for(j = 1; j <= ls; j++)
            {
                v = i * il + j;
                is = (i + j) * m + k;
                d[u] = d[u] - b[v] * d[is];
            }
        }
        if(ls != (il - 1)) { ls = ls + 1; }
    }
    return(2);
}

/**
 * @brief cgldl 求解对称方程组的分解法
 * @param a  n*n数组，存放方程组的系数矩阵（应为对称矩阵）
 * @param n  方程组的阶数
 * @param m  方程组右端常数向量的个数
 * @param c  n*m 存放方程组右端的m组常数向量；返回时存放方程组的m组解向量
 * @return   返回值小于0表示程序工作失败，大于0表示正常返回
 */
int cgldl(double a[], int n, int m, double c[])
{
    int i, j, l, k, u, v, w, k1, k2, k3;
    double p;
    if(fabs(a[0]) + 1.0 == 1.0)
    {
        return(-2);
    }
    for(i = 1; i <= n - 1; i++)
    {
        u = i * n;
        a[u] = a[u] / a[0];
    }
    for(i = 1; i <= n - 2; i++)
    {
        u = i * n + i;
        for(j = 1; j <= i; j++)
        {
            v = i * n + j - 1;
            l = (j - 1) * n + j - 1;
            a[u] = a[u] - a[v] * a[v] * a[l];
        }
        p = a[u];
        if(fabs(p) + 1.0 == 1.0)
        {
            return(-2);
        }
        for(k = i + 1; k <= n - 1; k++)
        {
            u = k * n + i;
            for(j = 1; j <= i; j++)
            {
                v = k * n + j - 1;
                l = i * n + j - 1;
                w = (j - 1) * n + j - 1;
                a[u] = a[u] - a[v] * a[l] * a[w];
            }
            a[u] = a[u] / p;
        }
    }
    u = n * n - 1;
    for(j = 1; j <= n - 1; j++)
    {
        v = (n - 1) * n + j - 1;
        w = (j - 1) * n + j - 1;
        a[u] = a[u] - a[v] * a[v] * a[w];
    }
    p = a[u];
    if(fabs(p) + 1.0 == 1.0)
    {
        return(-2);
    }
    for(j = 0; j <= m - 1; j++)
    {
        for(i = 1; i <= n - 1; i++)
        {
            u = i * m + j;
            for(k = 1; k <= i; k++)
            {
                v = i * n + k - 1;
                w = (k - 1) * m + j;
                c[u] = c[u] - a[v] * c[w];
            }
        }
    }
    for(i = 1; i <= n - 1; i++)
    {
        u = (i - 1) * n + i - 1;
        for(j = i; j <= n - 1; j++)
        {
            v = (i - 1) * n + j;
            w = j * n + i - 1;
            a[v] = a[u] * a[w];
        }
    }
    for(j = 0; j <= m - 1; j++)
    {
        u = (n - 1) * m + j;
        c[u] = c[u] / p;
        for(k = 1; k <= n - 1; k++)
        {
            k1 = n - k;
            k3 = k1 - 1;
            u = k3 * m + j;
            for(k2 = k1; k2 <= n - 1; k2++)
            {
                v = k3 * n + k2;
                w = k2 * m + j;
                c[u] = c[u] - a[v] * c[w];
            }
            c[u] = c[u] / a[k3 * n + k3];
        }
    }
    return(2);
}

/**
 * @brief chchol   求解对称正定方程组的平方根法, 用 Cholesky 分解法（即平方根法）
 * @param a n*n 存放系数矩阵（应为对称正定矩阵），返回时其上三角部分存放分解后的U矩阵
 * @param n 方程的阶数
 * @param m 方程组右端的常数向量的个数
 * @param d n*m 存放方程组右端m组常数向量；返回时，存放方程组的m组解向量
 * @return  返回值小于0，表示程序工作失败
 */
int chchol(double a[], int n, int m, double d[])
{
    int i, j, k, u, v;
    if((a[0] + 1.0 == 1.0) || (a[0] < 0.0))
    {
        return(-2);
    }
    a[0] = sqrt(a[0]);
    for(j = 1; j <= n - 1; j++)
    {
        a[j] = a[j] / a[0];
    }
    for(i = 1; i <= n - 1; i++)
    {
        u = i * n + i;
        for(j = 1; j <= i; j++)
        {
            v = (j - 1) * n + i;
            a[u] = a[u] - a[v] * a[v];
        }
        if((a[u] + 1.0 == 1.0) || (a[u] < 0.0))
        {
            return(-2);
        }
        a[u] = sqrt(a[u]);
        if(i != (n - 1))
        {
            for(j = i + 1; j <= n - 1; j++)
            {
                v = i * n + j;
                for(k = 1; k <= i; k++)
                {
                    a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j];
                }
                a[v] = a[v] / a[u];
            }
        }
    }
    for(j = 0; j <= m - 1; j++)
    {
        d[j] = d[j] / a[0];
        for(i = 1; i <= n - 1; i++)
        {
            u = i * n + i;
            v = i * m + j;
            for(k = 1; k <= i; k++)
            {
                d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j];
            }
            d[v] = d[v] / a[u];
        }
    }
    for(j = 0; j <= m - 1; j++)
    {
        u = (n - 1) * m + j;
        d[u] = d[u] / a[n * n - 1];
        for(k = n - 1; k >= 1; k--)
        {
            u = (k - 1) * m + j;
            for(i = k; i <= n - 1; i++)
            {
                v = (k - 1) * n + i;
                d[u] = d[u] - a[v] * d[i * m + j];
            }
            v = (k - 1) * n + k - 1;
            d[u] = d[u] / a[v];
        }
    }
    return(2);
}

/**
 * @brief ciggj  求解大型稀疏方程组的全选主元高斯-约当消去法,
 *               只是在消去过程中避免了对零元素的运算，从而可以大大减少运算次数.
 *               本函数不考虑系数矩阵的压缩存储问题
 * @param a n*n数组，存放方程组的系数矩阵，返回时要被破坏
 * @param n 方程组的阶数
 * @param b 存放方程组右端常数向量，返回时存放解向量
 * @param perm
 * @return 返回值为0表示系数矩阵奇异
 */
int ciggj(double a[], int n, double b[], int perm[])
{
    int i, j, k, is, u, v;
    double d, t;
    for(k = 0; k <= n - 1; k++)
    {
        d = 0.0;
        for(i = k; i <= n - 1; i++)
        {
            for(j = k; j <= n - 1; j++)
            {
                t = fabs(a[i * n + j]);
                if(t > d)
                {
                    d = t;
                    perm[k] = j;
                    is = i;
                }
            }
        }
        if(d + 1.0 == 1.0)
        {
            return(0);
        }
        if(is != k)
        {
            for(j = k; j <= n - 1; j++)
            {
                u = k * n + j;
                v = is * n + j;
                t = a[u];
                a[u] = a[v];
                a[v] = t;
            }
            t = b[k];
            b[k] = b[is];
            b[is] = t;
        }
        if(perm[k] != k)
        {
            for(i = 0; i <= n - 1; i++)
            {
                u = i * n + k;
                v = i * n + perm[k];
                t = a[u];
                a[u] = a[v];
                a[v] = t;
            }
        }
        t = a[k * n + k];
        for(j = k + 1; j <= n - 1; j++)
        {
            u = k * n + j;
            if(a[u] != 0.0)
            {
                a[u] = a[u] / t;
            }
        }
        b[k] = b[k] / t;
        for(j = k + 1; j <= n - 1; j++)
        {
            u = k * n + j;
            if(a[u] != 0.0)
            {
                for(i = 0; i <= n - 1; i++)
                {
                    v = i * n + k;
                    if((i != k) && (a[v] != 0.0))
                    {
                        is = i * n + j;
                        a[is] = a[is] - a[v] * a[u];
                    }
                }
            }
        }
        for(i = 0; i <= n - 1; i++)
        {
            u = i * n + k;
            if((i != k) && (a[u] != 0.0))
            {
                b[i] = b[i] - a[u] * b[k];
            }
        }
    }
    for(k = n - 1; k >= 0; k--)
    {
        if(k != perm[k])
        {
            t = b[k];
            b[k] = b[perm[k]];
            b[perm[k]] = t;
        }
    }
    return(1);
}

/**
 * @brief cjtlvs  求解 Toeplitz 型方程组的 Levinson 递推算法
 * @param t  长度为 n 的一维数组，存放对称T型矩阵中的元素t0, t1, ... tn-1
 * @param n  方程组的阶数
 * @param b  长度为 n 的一维数组, 存放方程组右端的常数向量
 * @param x  长度为 n 的一维数组, 返回方程组的解
 * @param s  长度为 n 的一维数组, 辅助用
 * @param y  长度为 n 的一维数组, 辅助用
 * @return  返回值小于 0 表示程序工作失败
 */
int cjtlvs(double t[], int n, double b[], double x[], double s[], double y[])
{
    int i, j, k;
    double a, beta, q, c, h;
    a = t[0];
    if(fabs(a) + 1.0 == 1.0)
    {
        return(-1);
    }
    y[0] = 1.0;
    x[0] = b[0] / a;
    for(k = 1; k <= n - 1; k++)
    {
        beta = 0.0;
        q = 0.0;
        for(j = 0; j <= k - 1; j++)
        {
            beta = beta + y[j] * t[j + 1];
            q = q + x[j] * t[k - j];
        }
        if(fabs(a) + 1.0 == 1.0)
        {
            return(-1);
        }
        c = -beta / a;
        s[0] = c * y[k - 1];
        y[k] = y[k - 1];
        if(k != 1)
        {
            for(i = 1; i <= k - 1; i++)
            {
                s[i] = y[i - 1] + c * y[k - i - 1];
            }
        }
        a = a + c * beta;
        if(fabs(a) + 1.0 == 1.0)
        {
            return(-1);
        }
        h = (b[k] - q) / a;
        for(i = 0; i <= k - 1; i++)
        {
            x[i] = x[i] + h * s[i];
            y[i] = s[i];
        }
        x[k] = h * y[k];
    }
    return(1);
}
